Phonon instability and self-organized structures in multi-layer stacks of confined 
dipolar Bose-Einstein condensates in optical lattices 
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In calculations to date [l|, [2| of multi-layer stacks of dipolar condensates, created in one- 
dimensional optical lattices, the condensates have been assumed to be two-dimensional. In a real 
experiment, however, the condensates do not extend to infinity in the oblate direction, but have 
to be confined by a trap potential, too. By three-dimensional numerical simulations of this realis- 
tic experimental situation we find a crucial dependence of the phonon instability boundary on the 
number of layers. Moreover, near the boundary of the phonon instability, a variety of structured 
ground-state wave functions emerges, which may indicate the onset of a roton instability [1, 

PACS numbers: 03.75.Lm, 03.75.Hh, 67.85.-d 



INTRODUCTION 

The experimental realization of Bose-Einstein conden- 
sates of chromium atoms [3], with a strong interatomic 
magnetic dipole-dipole interaction, has opened the door 
to promising new experiments Q on dipolar gases in 
which a wealth of new phenomena predicted by theory, 
such as the appearance of radial and angular rotons, bi- 
concave shapes of the ground state, anisotropic solitons, 
etc. 
able. 
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JJ, IM lla, should be observ- 
In dipolar condensates, the manipulation of the 



s-wave scattering length via Feshbach resonances offers 
the opportunity of tuning the relative strength of the in- 
teractions through the whole range, from dominance of 
the short-range interaction to that of the dipole-dipole 
interaction, and of studying the condensates in the dif- 
ferent regimes. 

Recently, stacks of non-overlapping dipolar conden- 
sates, unconfined in two dimensions, in one-dimensional 
optical lattices have been investigated theoretically P, Q . 
Special attention has been payed to the well-known 
roton-instability, and it has been reported in Ref. [H that 
this instability is enhanced, in the sense that the conden- 
sates become dynamically unstable in a wider range of 
the s-wave scattering length (measured in units of the 
dipole length) when more layers are added to the stack. 
Faraday pattern size transitions have been proposed 
as an experimental probe for the onset of the roton min- 
imum in multi-layer stacks of two-dimensional dipolar 
condensates. The phonon instability, however, in those 
calculations always set in at the same value of the ratio 
of contact and dipole-dipole interaction, independent of 
the number of layers. 

In an experiment a confining trap potential inevitably 
has to be present also in the planes of the conden- 
sates. It is the purpose of this paper to carry out three- 
dimensional numerical simulations for this realistic sit- 
uation, investigating the phonon instability in more de- 
tail. It will turn out that the latter depends crucially 
on the stack size. A second consequence of the three- 
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FIG. 1: A stack of Ns dipolar Bose-Einstein condensates. 
The long-range nature of the dipole-dipole interaction cou- 
ples all individual condensates. The traps are assumed to be 
extremely oblate and inter-site hopping of particles is prohib- 
ited. 



dimensional conhnement is the appearance of a variety 
of structured ground states in the vicinity of the phonon 
instability. It has been suggested in the literature 0, 01 
that the emergence of structures in the condensate den- 
sity is related to the onset of a roton instability. 



THEORY 

As in Ref. Q we consider a stack of Ng dipolar Bose- 
Einstein condensates, each containing the same number 
of atoms N. We consider axisymmetric traps, whose 
geometry is specified [l7| by the mean trap frequency 
uj — (cjgCJz)"'^/'^ and the aspect ratio A — LUz/t-Og {i-^g,z 
are the trap frequencies in the g- and z-direction, respec- 
tively). The spacing between adjacent traps. A, must be 
large enough for condensate wave functions not to over- 
lap (i.e. inter-site hopping of atoms is forbidden), but 
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sufficiently small for neighbouring condensates to "feel" 
each other, despite the relative weakness of the dipolar 
interaction. Hence the traps are extremely oblate (aspect 
ratios A > 100). Fig. [l]shows a schematic representation 
of the system. 

As is well known, at sufficiently low temperatures the 
ground state of a single condensate of weakly interact- 
ing bosons can be represented by a single wave function, 
whose dynamics obeys the Gross-Pitaevskii (GP) equa- 
tion. For a stack of condensates the wave functions in 
different layers are described by a set of coupled extended 
GP equations. To obtain this set of equations in dimen- 



sionless form, we introduce "atomic" units derived 
from the mass m and the magnetic dipole moment fi of an 
atom, namely a "dipole length" — TO/ioM^/(27r/i^), a 
unit energy = h?/ (2ma^), a unit frequency uj^ = Ej^/h 
and a unit time = h/E^, respectively. In addition, to 
eliminate the explicit dependence of the interaction po- 
tentials on the number of atoms, we carry out a particle 
number scaling for the position and time coordinates, 
r ^0.(1^ 1 1 ^ ^^^d^- scaled atomic units, the 

coupled system of GP equations for the set of Ns con- 
densate wave functions {^j{r,t)} then reads: 
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Hoj + 8™ \i}j{r, t)\^ + I dr 

\ k=l '' 



1-3 cos2 



r — r' 



|Vfe(r',Or i^j{r,t)=tdt^j{r,t), (j = 1, . 



(1) 



The kinetic energy and the trap potential terms for the 
condensate in layer j have been subsumed in iJoj , 



oj 



4„,2 



l-2j 



(2) 

where 7^.^ = ujg^z/i'^'-^d) denote the dimensionless trap 
frequencies. It is evident from ((!]) and ^ that the three 
external parameters governing the system are the scat- 
tering length a, measured in units of the dipole length a^, 
and the atom number scaled trap frequencies N'^^q.z, or, 
alternatively, the aspect ratio A — ^qI^z and the atom 
number scaled mean trap frequency A^^7 = N'^^j^J'^ "fV^ ■ 
We note that our parameter N'^^g is related to the di- 
mensionless parameter D used in Refs. [3, to mea- 
sure the effective strength of the dipole interaction in 
trap units by £> = (N'^^gj '^'^^'^ . We also note that the 
dimensionless parameter (i used by Klawunn and Santos 
Q to characterize the strength of the dipole interaction 



compared to the short-range interaction is related to our 
dimensionless scattering length a by /3 = l/(3a). 

Following Ref. Q , we make use of the very large aspect 
ratio A which leads to a decoupling of the motions in the 
longitudinal and the transverse direction 



(3) 



and assume that in the longitudinal direction the con- 
densate wave function is frozen in the ground state of 
the harmonic z-trapping potential in layer j 

0,(z,t) = ^-i/4;-i/2g-(.+^^-^A)V2P.g-.M.t^ (4) 

where fiz = N'^lz is the chemical potential of the ground 
state and Iz = i/l/Tt^ the oscillator length. Multiplying 
the GP equation ^ by (f>j and integrating over the z- 
coordinate yields after a lengthy calculation 
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o + Y^ \i'±A(^^t)\' {ji^{\lb±A9,t)f}Fikg,j-k)} J ib^Af^^t) = ^^ti^±A(^^t)^ (5) 



k=l 



with 



and 



F{kg, 1^) = ^ I dkz 



Ho = -A, + N^y, 



(6) 
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Hg -h K,z 



1 ) e-^"''^ cos(fc^nA). 

(7) 
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FIG. 2: Critical scattering length for the onset of roton- and 
phonon-instability depending on the number of layers. The 
mean trap frequency is = 1.3 ■ 10^. 
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FIG. 3: Critical scattering length for the onset of roton- and 
phonon-instability depending on the number of layers. The 
mean trap frequency is N^^ = 1.3 • 10*. 



In (l5|), and denote the Hankel transform and 

its inverse, respectively, which appear as a consequence 
of the 2D Fourier transform in cylinder coordinates used 
to evaluate the convolution integral. The ground state 
is obtained by solving ([5]) in imaginary time on a dis- 
cretized spatial grid. In the calculations we make use of 
the discrete Hankel transform method presented in [l9| . 

RESULTS AND DISCUSSION 

We simulate an experiment of a dipolar condensate of 
^^Cr atoms, for which ^ = and = 91a^, with 
and the Bohr radius and the Bohr magneton, respec- 
tively. Fig. [2] shows the results for the stability diagram 
for a mean trap frequency 7V^7 = 1.3-10^, a trap as- 
pect ratio A = 340 and a spacing A = 0.035. Choosing 
uj = 27r • 1.0 kHz, this implies approximately 3300 parti- 
cles per layer and a spacing of 0.56 |J.m. These are typical 
numbers for actual experiments. It is evident that there 
appears a strong dependence of the critical value of (3 
(or, in other words, the scattering length a, with a < 0) 
on the number of layers. This instability threshold is 
denoted phonon instability due to the fact that excita- 
tions of the condensates with small or vanishing wave vec- 
tor (phonons) lead to imaginary energy eigenvalues and 
consequently collapse of the condensates. The fact that 
the phonon instability moves to larger absolute values of 
(3 (i.e. smaller absolute values of the scattering length) 
when the number of condensates is increased can be un- 
derstood in a very simple picture: With growing number 
of sites, the total dipolar potential to which the indi- 
vidual condensates are exposed becomes stronger. This 
potential forces the condensates to shrink radially which 
leads to a higher density, necessary for inducing collapse. 
Consequently, the collapse threshold is being pushed to 
larger absolute values of (3. As our calculations show, it is 
the condensate in the middle of the stack which becomes 
unstable first. This is possibly not surprising, since the 
dipolar potential approaches a maximum in the middle 



layer. The phonon instability must therefore be under- 
stood in the sense that below the threshold, there is at 
least one layer, and not necessarily the whole system, 
which becomes unstable. 

At Ns ~ 50, /?crit approaches a constant value. Again, 
we can explain this by considering the condensate in the 
middle. The dipole-dipole interaction decays as l/r'^, and 
the sum over the dipolar potentials of all condensates 
converges with growing number of sites, as well as the 
phonon instability. We point out that the effect described 
can only appear in a 3D trap. In 2], the homogeneous 
case was considered, where the condensate density was 
the same for all layers and for any number of sites. As 
a consequence, the phonon instability appeared at the 
same value of (3 independent of the number of sites. 

Fig. [3] shows the results for the stability diagram for 
a mean trap frequency of N"^^ — 1.3-10^ and a value of 
the spacing parameter of A = 0.011. The latter value 
is reduced by a factor of a/TO because we effectively in- 
creased the particle number in the system by a factor of 
a/IO. According to the scaling laws, in the experiment 
proposed, the spacing is still 0.56 (a.m. As we can see in 
Fig. [HI in the limit Ng ^ 00, (3ait approaches a larger 
value than in Fig. [2l This is again due to the higher 
density in the condensates. 

Another feature that already appears in single con- 
densates is the formation of structures in the conden- 
sate wave function. It has been reported [1, |3| that near 
the phonon instability and at special trap aspect ratios, 
the shape of the condensate wave function differs from 
a "normal" Gaussian or parabola, where the maximum 
density lies at the centre of the trap, and instead forms 
a density pattern with the peak density away from the 
centre. Due to the special shape these density patterns 
have been dubbed as "blood-cell-like" condensates. 

In extremely oblate and layered condensates, the situ- 
ation is qualitatively different. Fig. shows the radial 
part of the wave function for a single condensate trapped 
at A^^7 = 1.3 • 10^ in the very vicinity of the phonon 
instability. One peak of the density distribution lies at 
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FIG. 4: Structures of the wave functions in the very vicinity 
of the phonon instability for a single trap (a) and 2 traps (b). 
The mean trap frequency is — 1.3 • 10^. 
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FIG. 5: Structures of the wave functions in the very vicinity 
of the phonon instability for 4 traps (a) and 5 traps (b). The 
mean trap frequency is N^y = 1.3 • 10^. 

the centre of the trap, but there is also a second, larger 
peak at g « 0.09. This picture changes when a second 
layer is added to the stack. The radial wave functions 
for the condensates 1 and 2 (which are identical due to 
the symmetry) now assume the familiar blood-cell shape 
with one density peak away from the centre of the trap. 



This structure remains when further layers are added. In 
the stability diagrams, Figs. [2] and [Sj the regions where 
these structures arise are coloured blue. When increasing 
the trapping frequency to N'^^ — 1.3 • 10*, the situation is 
even richer. Fig. [5^ shows the radial wave functions for a 
stack of 4 condensates. In all layers, the condensate wave 
functions exhibit 2 peaks away from the centre, while the 
peaks are more pronounced for the two in-most conden- 
sates; this shows again the effect of the dipole-dipole in- 
teraction being enhanced by the outer condensates. In 
Fig. [SJa, a fifth layer was added, and it now seems en- 
ergetically more favourable for the condensates to move 
one peak to the centre, while keeping an outer, smaller 
peak. 



CONCLUSIONS 

By three-dimensional numerical simulations of multi- 
layer stacks of confined dipolar Bose-Einstein conden- 
sates in optical lattices for realistic experimental parame- 
ters we have revealed a crucial dependence of the phonon 
instability boundary on the number of layers. This find- 
ing is at variance with the results of Ref. 3 where un- 
confined 2D condensates were considered and the criti- 
cal value of the boundary was independent of the num- 
ber of layers. We have also demonstrated that near the 
boundary of the phonon instability, a variety of struc- 
tured ground-state wave functions emerges. 

Structured ground-state wave functions have been in- 
tensively discussed in the literature [1, [3l- It has been 
shown that the emergence of structures is intimately re- 
lated to the appearance of imaginary energy eigenvalues 
in the Bogoliubov spectrum. More precisely, the Bogoli- 
ubov modes that become unstable have non-vanishing 
projection of the angular momentum on the 2;-axis [ij. 
These angular roton modes are responsible for a dynam- 
ical instability of the condensates resulting in local col- 
lapse. Contrary to global collapse, where the particles are 
accumulated at the centre of the trap, a locally collaps- 
ing condensate forms several density peaks ordered on a 
ring. It is the same ring on which the ground-state wave 
function reaches a maximum. In our case, the situation 
is different since more than one density peak occurs. Ev- 
idently a 3D numerical investigation of the Bogoliubov 
spectrum along these lines in the region of structured 
ground states in Figs. [5] and [3] is urgent, and under way. 

We thank Jonas Metz and Stefan Miiller from 5. 
Physikalisches Institut, Universitat Stuttgart, for fruit- 
ful discussions. 
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